## Copyright (C) 2001 Paulo Neis
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, write to the Free Software
## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
## 

## N-ellip 0.2.1
##usage: [Zz, Zp, Zg] = ellip(n, Rp, Rs, Wp, stype,'s')
##
## Generate an Elliptic or Cauer filter (discrete and contnuious).
## 
## [b,a] = ellip(n, Rp, Rs, Wp)
##  low pass filter with order n, cutoff pi*Wp radians, Rp decibels 
##  of ripple in the passband and a stopband Rs decibels down.
##
## [b,a] = ellip(n, Rp, Rs, Wp, 'high')
##  high pass filter with cutoff pi*Wp...
##
## [b,a] = ellip(n, Rp, Rs, [Wl, Wh])
##  band pass filter with band pass edges pi*Wl and pi*Wh ...
##
## [b,a] = ellip(n, Rp, Rs, [Wl, Wh], 'stop')
##  band reject filter with edges pi*Wl and pi*Wh, ...
##
## [z,p,g] = ellip(...)
##  return filter as zero-pole-gain.
##
## [...] = ellip(...,'s')
##     return a Laplace space filter, W can be larger than 1.
## 
## [a,b,c,d] = ellip(...)
##  return  state-space matrices 
##
## References: 
##
## - Oppenheim, Alan V., Discrete Time Signal Processing, Hardcover, 1999.
## - Parente Ribeiro, E., Notas de aula da disciplina TE498 -  Processamento 
##   Digital de Sinais, UFPR, 2001/2002.
## - Kienzle, Paul, functions from Octave-Forge, 1999 (http://octave.sf.net).
##
## Author: Paulo Neis <p_neis@yahoo.com.br>
## Modified: Doug Stewart Feb. 2003


ellip <- function(n, ...) UseMethod("ellip")

ellip.FilterOfOrder <- function(n, Rp = n$Rp, Rs = n$Rs, W = n$Wc, type = n$type, ...)
  ellip(n$n, Rp, Rs, W, type, ...)

ellip.default <- function(n, Rp, Rs, W, type = c("low", "high", "stop", "pass"), plane = c("z", "s"), ...) {

  type = match.arg(type)
  plane = match.arg(plane)

  ## interpret the input parameters
  if (!(length(n)==1 && n == round(n) && n > 0))
    stop("ellip: filter order n must be a positive integer")

  stop = type == "stop" || type == "high"
  digital = plane == "z"

  if (length(W) != 1 && length(W) != 2)
    stop("ellip: frequency must be given as w0 or c(w0, w1)")

  if (digital && !all(W >= 0 & W <= 1))
    stop("ellip: critical frequencies must be in (0 1)")
  else if (!digital && !all(W >= 0))
    stop("ellip: critical frequencies must be in (0 inf)")

  if (Rp < 0)
    stop("ellip: passband ripple must be positive decibels")

  if (Rs < 0)
    stop("ellip: stopband ripple must be positive decibels")

  
  ## Prewarp to the band edges to s plane
  if (digital) {
    T = 2       # sampling frequency of 2 Hz
    W = 2 / T*tan(pi * W / T)
  }

  ##Generate s-plane poles, zeros and gain
  ZPG = ncauer(Rp, Rs, n)

  ## s-plane frequency transform
  ZPG = sftrans(ZPG, W = W, stop = stop)

  ## Use bilinear transform to convert poles to the z plane
  if (digital)
     ZPG = bilinear(ZPG, T = T)

  as.Arma(ZPG)
}

ncauer <- function(Rp, Rs, n)  { #  [zer, pol, T0]

  ellip_ws <- function(n, rp, rs) {
    ellip_ws_min <- function(kl) {
      int = ellipke(c(kl, 1-kl))$k
      ql = int[1]
      q = int[2]
      abs((ql/q) - x)
    }
    kl0 = ((10^(0.1*rp) - 1) / (10^(0.1*rs) - 1))
    k0 = 1 - kl0
    int = ellipke(c(kl0, k0))$k
    ql0 = int[1]
    q0 = int[2]
    x = n * ql0 / q0
#    kl = optim(ellip_ws_min, method = "L-BFGS-B", lower = eps, upper = 1-eps)
    kl = optimize(ellip_ws_min, interval = c(eps, 1-eps))$minimum
    ws = sqrt(1/kl)
    ws
  }

  ## Cutoff frequency = 1:
  wp = 1
  
  ## Stop band edge ws:
  ws = ellip_ws(n, Rp, Rs)
  
  k  = wp / ws
  k1 = sqrt(1-k^2)
  q0 = (1/2)*((1-sqrt(k1))/(1+sqrt(k1)))
  q =  q0 + 2*q0^5 + 15*q0^9 + 150*q0^13

  ####Filter order maybe this, but not used now:
  ##D =  (10^(0.1*Rs)-1)/(10^(0.1*Rp)-1)
  ##n=ceil(log10(16*D)/log10(1/q))

  l=(1/(2*n))*log((10^(0.05*Rp)+1)/(10^(0.05*Rp)-1))
  sig01=0
  sig02=0
  for ( m in 0:30) {
	sig01=sig01+(-1)^m * q^(m*(m+1)) * sinh((2*m+1)*l)
  }
  for ( m in 1:30) {
	sig02=sig02+(-1)^m * q^(m^2) * cosh(2*m*l)
  } #
  sig0=abs((2*q^(1/4)*sig01)/(1+2*sig02))

  w=sqrt((1+k*sig0^2)*(1+sig0^2/k))
  
  if (n %% 2) # odd
	r=(n-1)/2
  else
	r=n/2
  
  wi=matrix(0, 1,r)
  for (ii in 1:r) {
	if (n %% 2) #odd
      mu=ii
	else
      mu=ii-1/2
	soma1=0
	for (m in 0:30) {
      soma1 = soma1 + 2*q^(1/4) * ((-1)^m * q^(m*(m+1)) * sin(((2*m+1)*pi*mu)/n))
	}
	soma2=0
	for (m in 1:30) {
      soma2 = soma2 + 2*((-1)^m * q^(m^2) * cos((2*m*pi*mu)/n))
	}
	wi[ii] = soma1 / (1+soma2)
  }
  
  Vi = sqrt((1-(k*(wi^2)))*(1-(wi^2)/k))
  A0i = 1 / wi^2
  sqrA0i = 1 / wi
  B0i = ((sig0*Vi)^2 + (w*wi)^2) / ((1+sig0^2*wi^2)^2)
  ## not used: ## B1i = (2 * sig0 * Vi) / (1 + sig0^2 * wi^2)

  ##Gain T0:
  if (n %% 2) # odd
	T0=sig0 * prod(B0i / A0i) * sqrt(ws)
  else
	T0=10^(-0.05*Rp) * prod(B0i / A0i)
  
  ##zeros:
  zer = c(1i*sqrA0i, -1i*sqrA0i)

  ##poles:
  pol =c((-2*sig0*Vi + 2*1i*wi*w) / (2*(1 + sig0^2*wi^2)), (-2*sig0*Vi - 2*1i*wi*w) / (2*(1 + sig0^2*wi^2)))
  
  ##If n odd, there is a real pole  -sig0:
  if (n %% 2) # odd
    pol = c(pol, -sig0)

  ##
  pole = sqrt(ws) * pol
  zero = sqrt(ws) * zer

  Zpg(zero = zero, pole = pole, gain = T0)
} 
